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Abstract. In this article we describe applications of the numerical method of 
discrete differential forms in computational GR. In particular we consider the 
initial value problem for vacuum space-times that admit plane gravitational 
waves. As described in an earlier paper the discrete differential form approach 
can provide accurate results in spherically symmetric space-times [28]. More- 
over it is manifestly coordinate independent. 

Here we use the polarised Gowdy solution as a testbed for two numerical 
schemes. One scheme reproduces that solution very well, in particular it is 
stable for a comparatively long time and converges quadratically. 



1. Introduction 

In an earlier paper on the application of discrete differential forms in numeri- 
cal General Relativity we described this method for spherically symmetric prob- 
lems [28]. Here we apply the numerical schemes in the context of symmetric 
cosniological space-times. 

In [28] we used discrete differential forms to calculate the geometry of the do- 
main of dependence of a given initial hypersurface. In that context the domain of 
dependence was only a small section of the full space-time, which makes it impos- 
sible to investigate the longterm behaviour of the numerical schemes unless one is 
prepared to introduce outgoing boundary conditions. However, this would not al- 
low us to judge the pure evolution algorithm. Thus, we consider space-times which 
have compact Cauchy surfaces, i.e. cosmological solutions of Einstein's equations. 

Thus, we will consider space-times M with the topology S x R, where S is a 
closed 3-manifold homeomorphic to the Cauchy surfaces of Ai. The topology of 
the Cauchy surface can be controlled quite easily within the discrete differential 
form approach, because it is given in the very beginning in terms of incidence 
matrices [6] . Hence, we do not need to worry about outgoing boundary conditions 
anymore. For instance, if the topology of S is the 3-torus then one may equally 
well use the interpretation that periodic boundary conditions are imposed in each 
direction. We will not consider the full 3-dimensional case here but as in [28] we 
will confine ourselves to cases where the problem is effectively 1 -I- 1-dimensional by 
imposing appropriate symmetries. 

An example for a symmetric space-time which admits Cauchy surfaces with 
topology T'^ is the polarised Gowdy solution [15] . This solution can be interpreted as 
describing gravitational waves in a cosmological universe. In contrast to the static 
spherically symmetric space-times of [28] this solution is time-dependent. Using 
this solution to test the method has another advantage: it is one of the testbeds 
for numerical codes suggested by the Apples with Apples alliance [1] and thus one 
may compare the results with other numerical methods. 
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The plan of the article is as follows. In section [2] we shortly summarise the 
properties of discrete differential forms (more information can be found in [6, 10, 
28]). Then, in section [3] we describe the equations which result from a symmetry 
reduction. In section [4] we present the method of implementing these equations in 
two fully discrete numerical schemes. In section[5]we discuss how the schemes were 
tested, and in section [6] we present the results of those tests. 

2. Preliminaries 

Discrete differential forms in numerical GR were first introduced in [10]. The 
main feature of this method is its manifest coordinate independence. Coordinate 
invariance is one of the fundamental properties of GR. Thus, it is natural to use 
coordinate invariant numerical methods. A discretisation procedure with that prop- 
erty is Regge calculus [13,27]. However, so far it did not play a role in numerical 
GR. In computational approaches to treat the problem of coordinate dependencies 
multiple coordinate systems are used to cover the space-time [14,30,31,35], but 
there the coordinate invariance is not manifest. 

A fundamental difference of the discrete differential form approach in compar- 
ison with many other numerical methods is the discretisation of space-time. In 
contrast to the usual procedure, where the finite counterpart of the manifold is a 
numerical grid that is composed of a finite number of points (see e.g. [2] for a review 
article about numerical GR), the discrete space-time here also contains objects of 
higher dimensionality like curves and surfaces. The collection of points (nodes), 
curves (edges), surfaces (faces) and volumes is a cellular paving [6] and it is called 
computational mesh. The various elements of a cellular paving are called cells and 
through incidence relations between the cells the topology of the computational do- 
main is defined. We use cellular pavings whose elements are all simplices, because 
this provides an elegant way to define a discrete version of the exterior product 
(see [10,28]). The cellular paving then corresponds to a simplicial complex [19]. 

Based on the cellular paving one discretises a theory that is formulated in terms 
of differential forms. Differential p-forms can be viewed as 'the objects which are 
integrated over p-dimensional submanifolds'. That means they provide maps from 
p-dimensional submanifolds to the reals. Then discrete p-forms are maps that 
assign a number to every p-dimensional cell in the cellular paving. They have 
received some attention since Bossavit [5] had pointed out that they correspond to 
the lowest order mixed finite element spaces defined by Nedelec [20] (see also [26]). 
Finite elements of mixed type have been used successfully in numerical applications 
to electrodynamics, see [4,6,16]. In numerical GR the finite element method has 
been used just recently e.g. in [32]. 

In order to use this approach one needs to have a formulation of geometries and, 
in particular, of GR which uses differential forms. A formulation of geometries based 
on differential forms has been provided by E. Cartan [7]. The further step towards a 
formulation of GR using differential forms has been carried out by several authors. 
We mention here the work of Sparling [33] who has set up an exterior differential 
system of equations which is closed if and only if the vacuum Einstein equations 
hold. In [10] it is shown in detail how to set up the discrete formalism based on 
this exterior system using the ideas explained above. 

In summary, the variables of our proposed discrete formulation will be the inte- 
grals of the differential forms in Sparling's formulation of GR over submanifolds. 
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In order to get a finite number of variables we use a finite number of these sub- 
manifolds based on a triangulation of the computational domain. In this article we 
describe a simplification of the general formalism which occurs in space-times with 
a two-dimensional translational symmetry (see appendix E]) . 

3. The effective equations 

To obtain the effective equations in the space-times of interest we start with the 
Cartan formulation of GR (using exterior forms) [33] . The basic variables in this 
formalism are the four 1-forms of a pseudo-orthonormal tetrad 9^ , i — 0, ... ,3 [18]. 
With them the metric is 

(1) g = e^ -9^ -e"^ (gie"^ -6^ ^ rj.kO' ® 6/^ 

For the description of the connection in this formalism sixteen 1-forms u}^k, i,k — 
0, . . . ,3 are used. The connection should be compatible with the metric and tor- 
sion free, which translates to the antisymmetry-requirement and the first Cartan- 
equation respectivelj0: 

(2) 7/,fca;^ + r^jfcu;^ = 0, dO' + uj\e'' = 0. 

Furthermore, the metric should fulfil Einstein's vacuum field equations, which is 
equivalent to 

(3) E, - 0, 
where E^ is the Einstein 3- form defined by (see [33]) 

(4) E, = i£„fcif2^'"0\ 

with the curvature 2- form Cl^ k = du)^ k + (j^-'i^^^'k [9]. 

Since the antisymmetry of the connection 1-forms can easily be imposed, we are 
thus interested in the following system of equations 

(5a) d9' + uj\e^ ^ 0, 

(5b) E, = 0. 

Sparling considered these equations on the frame bundle over the space-time man- 
ifold. He showed that Einstein's equations are satisfied if and only if the ideal 
generated by ([51) is a closed differential ideal [33] . 

Even though the geometry is fixed by ([5]) , there is still the freedom of choosing a 
gauge, i.e. there are Lorentz transformations K^k of the tetrad that do not change 
the metric 

(6) g = ri^te' (8) = v^k{A',e^) ® (A^S') = {^kA^jA^) 6' ® oK 

In this work we will concentrate on general relativistic systems that occur in the 
context of space-times with a two-dimensional translational symmetry. Essentially 
that means that the group of translations (or T^) acts isometrically on the 
space-time and that the orbits of this action are two-dimensional, space-like, flat 
submanifolds (the 'wave fronts'). The details are discussed in appendix|3 where it 
is also shown how to 'factor out' the symmetry action and how to derive an exterior 
system on the two-dimensional space M.i oi orbits. 



^Here and in what follows it is understood, that the product of differential forms is the anti- 
symmetrised tensor-product, i.e. the exterior product. 
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This procedure yields two exterior systems. The first one is 

d(/i + 9i) - (/o + .90)'^ + (/o/i + 50.91)^° 
(7a) +(/i'-/o.9o + /igi+.9?)^' =0, 

d(/o + 50) - (/i + 91)^ + Uoh + go9i)0' 

(7b) +(/2 - A.91 + /ogo + .9o)^" = 0, 

(7c) da; + d(.go0' + 51^°) + (.9o - 9l)0°e' = 

(7d) da; + difo9' + /i©") + {f^ - fDO^O^ = 0, 

(7e) d6»" + a;6/i = 0, 

(7f) dO^ + ljG^ = 0, 

(7g) d (/o0° + h0^) - 0. 

Here (6>°,6>^) is a dyad in the two-dimensional orbit space Aii which carries a 
Lorentzian metric. The 5'0(1,1) connection on this space is given by the 1-form 
(jj'^i =: u3. It is a consequence of the equations above that this connection is torsion 
free. The geometric properties of the orbits are described by the functions /o, /i, 
50 and gi. For details see appendix VK\ 

The other exterior system that we shall consider is obtained by a procedure 
analogous to the spherically symmetric case [28] . For the system ([7]) it is performed 
in [29]. One ends up with the following system 

(8a) 

(8b) 

(8c) 

(8d) 

(8e) 

(8f) 

(8g) 

where the 1-forms a, (3, 7 and d are defined as 

a:=/o0" + /i0\ (3 := he" + he\ 

(9) 7:=go0"+5i0\ J := .gi^" + go©'. 
To impose these correlation one uses the following algebraic equations 

(10) = /3, *7 = 5, 

where ★ is the two-dimensional Hodge-operator [9] . 

The systems (0 and ([St.lfTU]) respectively now serve as the bases to derive nu- 
merical schemes. To achieve that we use the methods described in [6, 10,28]. 

4. Implementation of the discrete equations 

In section [3] and appendix |^ we derived two systems of equations on the orbit 
space M.I. Now we explain how these systems are discretised and develop numerical 
schemes. The first scheme is based on the system ([7]). It can be seen as the analogue 
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Figure 1: The triangulation S of the first four time-steps in the evolution of Ci. 



of scheme I in [28]. The second scheme is based on ((8t. (|T0)) using the same ideas as 
scheme III in [28]. 

4.1. Properties of the simplicial mesh. The first step in the discretisation 
procedure is the definition of the finite counterpart of the manifold. In our case 
this step is common to both schemes. As we mentioned in section [2l that structure 
is a simplicial complex S, and since we are in two dimensions that means it is 
composed of simplices of dimensions zero (nodes), one (edges) and two (faces). 

Here the simplicial approximation 5 of a subset of the orbit space A^i is gener- 
ated in analogy to the procedure described in [28]. Thus, the starting point for the 
construction of S is an initial hypersurface I and its discretisation, a 1-dimensional 
simplicial complex Ci [19] , which we define so that it approximates a space- like curve 
in A^i (see section [5^ for details). The difference here is that the topology of 2 is 
S^, i.e. that we impose periodic boundary conditions. The domain of dependence 
of the initial hypersurface is thus the full orbit space Aii- 

Having obtained Ci, a very elegant and invariant way to define the positions of 
the nodes in 5 \ Ci is to send lightrays h from the nodes in Ci and to take the 
intersections of these light-like geodesies as those positions. 

By identifying links between these nodes with the edges of the simplicial complex 
appropriately one obtains the simplicial approximation S. The construction is 
described in detail in [28] and its result is illustrated in figure [T] 

The complex S contains two types of faces, the upwards and the downwards di- 
rected ones. In figure[l]the downwards directed faces are hatched, and the upwards 
directed ones are not. Moreover it contains two types of edges, namely space-like 
and light-like ones. Collecting the space-like edges one obtains a structure that, on 
the level of topology, appears as N copies of the initial complex Ci (dashed lines in 
figure [1]) . Propagation from one of these copies (time-slice) to the next one is one 
time step. 

A p-dimensional simplex is commonly denoted [uq, . . . , n^], where the rii are the 
corners of the simplex. In what follows we will denote the upwards and downwards 
directed faces [no,ni,n2] and [rig, n'j^, n'j] respectively, such that the edge [^Iq'', ?t-i'''] 
is space- like and the edges [n^''', 71,2''] as well as [n-o\'T-2''] light-like (see figure 
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Figure 2: The upwards directed faces are denoted [no,ni,n2] and the downwards 
directed faces are [?1q, n'^, n^. 



As we already discussed in [28] the construction through Hght-like geodesies 
seems to be the simplest invariant method to define the position of 712 in (1 + 1)- 
dimensional manifolds. The choice of the nodes at a later time and their connections 
to the nodes at the initial time is essentially arbitrary and only restricted by topo- 
logical considerations. It is only when the 0' are known on all the edges that the 
geometry of the mesh is determined. We will see later that a part of these values 
can be specified freely while the rest is determined from the equations. 

As in all numerical simulations degeneracies may occur. For instance two ad- 
jacent nodes in the same time slice may have a time-like distance. However, this 
must be seen as a sign that the mesh is too coarse and should be refined. 



4.2. Properties of discrete forms. To discretise the equations ([7]) and ([5]) one 
needs discrete versions of the operations exterior product and exterior derivative. If 
the system pO|) is taken into account additionally a discrete Hodge operator must 
be defined. 

To get a discrete exterior product we use the method explained in [10] leading 
to the following formulas. Let a.^ and /3' be discrete p- and g-forms respectively 
then 

Q:°/3«[no, ...,%] = {a" [no] + ... + a°[ng]) /3«[no, . . . ,n,], 

(^11^ a^l3'^[no,ni,n2] = ^ (Q!^[no, ni]/3^ [no, 712] - a.^[no,n2](3^[no,ni] 

+ a^[ni, n2]/3^ [ni ,no]~ [ni , no]/3^ [ni , 712] 

+ a^[n2,no]f3^[n2,ni] - o;^ [n2, ni]/3^ [712, tiq]) . 

To obtain a natural discrete exterior derivative one applies Stokes' theorem. For 
us the relevant formulas are (see [6,10]) 

da"[no, rii] = Q!"[ni] — a"[no], 
(X2^ 1 111 

da [no, ni, n2] = a [ni,n2]-Q; [no,n2] + a [no,ni]. 

The Hodge operator is discretised as described in the context of scheme HI in [28] . 
This discretisation is based on the following observation. Let a and f3 be 1-forms 
such that a = -kfS. When this equation is integrated over a light-like geodesic 
with e° =e [e e {-1,+1}) then also 

(13) / a = e I 13. 

•lie J-fe 



DISCRETE DIFFERENTIAL FORMS FOR COSMOLOGICAL SPACE- TIMES 



7 



Thus wc use the foUowing discretisation of the Hodge operator. When the edge 
[ni,n2\ is hght-likc then 

(14) ★/3[ni,7i2] = £a[ni,n2]. 

In the Hodge operator is only apphed to l-forms. Therefore it is not necessary 
to define its discrete version for general forms here. 

Having a simplicial mesh and the discrete operators we can now develop the 
numerical schemes. Common to both schemes is that for each triangle a system 
of equations has to be solved. These are coupled non-linear algebraic equations. 
Their analysis is somewhat complicated and their properties are not yet clear. They 
might not have a unique solution. However, at least one solution can be found by 
Newton's iteration method. We used the GNU Scientific Library, especially the 
implementation of the modified Powell method [11,24]. 

4.3. Scheme I. The first scheme is based on the system ([7]). The variables are the 
discrete l-forms 0°, 0^ and u) as well as the discrete 0- forms /o, /i, go and gi. For 
the upwards directed faces the numbers 

{fo[no]Ja[ni], fi [tiq] , /i [ni ] , go [uq], go[ni], gi[no], gi[ni], 

(15) 0"[no, ni], 0^[no, ni], a;[no, ni]} 
are given initial data, and 

{/o ["-2] , /i [n-2] , go [«2] , 51 [^2] , 

(16) 0°[no,n2],0'^[ni,n2],0^[no,n2],6^[ni,n2],u}[nQ,n2],uj[ni,n2]} 
are the unknowns. 

The system ([7]) is composed of two 1-form equations d7a|) . (j7bp and five 2-form 
equations (I7cp-(|7g[). Therefore, since there are two 1-form equations on each of 
the light-like edges and five 2-form equations on the upwards directed face, it is 
clear that the total number of equations in the discretised system is nine. However, 
we discuss in appendix [X] that on the continuous manifold the 1-form equations 
([7a|) . (|7b[) take a special role, because they are not independent of the remaining 
equations in ([7]). They are partially redundant. Therefore we require the discrete 
1-form equations (|7a[) . (j7bp to be satisfied at the edge [no, ^2], but we do not require 
these equations at the edge [711,712]. 

Hence, we have seven equations for the upwards directed faces and ten unknowns. 
We get three more equations by using the definition of the position of 712 and fixing 
the gauge outside the initial hypersurface. That means we require 

(17) (0°-0l)[71o,7l2] -0, (0°+0')[71l,7l2] =0 

to specify the position of n2, and the gauge is chosen such that 

(18) u:[no,7i2] + <^[ni,n2] ^ 0. 

This gauge choice was also used in scheme HI of [28] . 
For the downwards directed faces the unknowns are 

(19) {0'[n'o,n[],0'[n'„n[],u;[n'o,n[]}. 

We already discussed how the observation, that the equations ([7a|) . (|7b|) are re- 
dundant, can be implemented on the discrete level. Following that reasoning we 
conclude that those 1-form equations need not be required at [no,n'i]. 
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It follows that we only need to consider the system ([7c|) -([7gl). These are five 
2-form equations, but the number of unknowns is only three. Thus, the system is 
overdetermined and we have to select three of the five equations (|7c|) -([7gl). At the 
moment it is not clear how a reasonable choice should be made. Compared to other 
choices that were tested shortly and provided rapidly growing errors the following 
system seems to work quite well 

(20a) = d6»" + u;6/\ 

(20b) O^de^+uj0°, 

(20c) = do; + d(5i0O + goO') + {gl - 5?)^°^'- 

Altogether these are seven equations and seven unknowns for the upwards di- 
rected faces, as well as three equations and three unknowns for the downwards 
directed ones. 



4.4. Scheme II. To obtain the second scheme we discretise the equations ([8al 
iroi), i.e. 



0^de°+u}e\ O = d0i+a;6»°, 

= dw + d^ + 7^, = da; + d/3 + af3, 

(21) = d7, = da, 
= d(^ + ^) + (a+7)(/3 + 5), 

= -k(3 — a, = ^(5 — 7. 

That means the variables are the seven 1-forms a, (3, 7, 9^ , 9^ and a;. The 
given initial data are hence 

{a [no, ni],/3[no, ni],7[no, rii], 5[no, "1], 

(22) 0°[no, ni], 0' S, rii], u;[no, ni]}, 
and the unknowns for the upwards directed faces are 

{a [no , 712] , a [ni , 712] , /3 [no , 712] , /3 [ni , 712] , 
7[no , n2] , 7 [ni , n2] , ^ [no , n2] , ^ [ni , 712] , 

(23) 0° [no , n2] , 0" [ni , n2] , 0^ [no ,712], 9^ [ni , n2] , a; [no , n2] , a; [ni , n2] } . 

The system (PT|) is composed of seven 2-form equations and two algebraic equa- 
tions involving the Hodge operator. These equations are discretised with the ex- 
terior product pT|) . the exterior derivative and the Hodge operator (|14p . For 
the upwards-directed faces this procedure results in seven equations for every face 
and two equations for both of the two light-like edges. In total these are eleven 
equations. 

With the same procedure as in scheme I, i.e. using that the new edges are light- 
like and choosing a gauge with a;[no,n2] -I- a;[ni,n2] = 0, we reduce the number 
of unknowns from fourteen to eleven, such that the number of equations for the 
upwards directed faces equals the number of unknowns there. 

For the downwards directed faces the unknowns are the integrals of the seven 
1-forms along the space-like edge, i.e. we have seven unknowns 

{a [no , n'l ] , /3 [nj, , n'l ] , 7 [n^ ,n[],5 [n^ , n'l] , 

(24) 9''[n'^,n[],9\n'^,n[],u:[n'^,n',]}. 
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The discrete equations involving the Hodge operator are aheady satisfied at the 
hght-hke edges [?^o,?^2] [71^,712]. Hence we get seven discrete equations from 
the 2-forms in ([H]), i.e. from (|5a))-([5g|). 



5. Test-scenarios 

In the last chapters we described, how discrete differential forms can be applied to 
study systems with a planar symmetry in GR. In this chapter we want to present 
the concrete example, where the code was tested. The idea for the test is to 
derive discrete initial data from an analytical solution, use the numerical schemes 
to simulate the time evolution and finally compare the numerical results with the 
analytically expected ones. That means an analytical solution is needed. 

As a testbed we have chosen the polarised Gowdy space-time [1, 15], because it 
has at least three advantages. Its spatial slices have the topology T'^, the light-like 
geodesies take a very simple form and it is suggested as a testbed for numerical 
schemes by the Apples with Apples alliance [1]. In standard coordinates {t, z, x, y} 
its metric reads 

(25) g = r^l'^e^l'^ {dt (g) di - dz dz) - te^dx (g) dx - te^^dy (g) dy, 
where i G R, z G [0, 1], 

P{t, z) ~ Jo{2TTt) cos27rz, 

X{t, z) = -27rt Jo(27rt) Ji(27rt) cos^ 2ttz + [^0 (Svrt) + Ji (27ri)] 

- 27r2 [Jo^(27r) + JI{2ti)\ + ttJq{2tt)Ji{2tt), 

and Jo, Ji are the usual Bessel- functions. 

It can easily be checked that a light-like geodesic in the {t, z)-surface through a 
point with coordinates {to, zq) also contains the points {to + k, zo) {k G N). In this 
sense, a massles test particle needs a time oi 5t = 1 to travel around the universe. 
This time interval is commonly denoted a crossing time. 



5.1. The continuous forms. To get the differential forms we make a gauge-choice, 
i.e. we choose some 0^ and 0^ , that generate the corresponding metric. In polarised 
Gowdy geometry a natural choice leads to 

/o = \t'^e--^ + dtP^ , A = \t--e--^d.P, 

go = ^t''e 4(--atPj, .91 = -2^*6 '^d^P, 

0° = t-^'^e^/^dt, 0' = t-^/^e^'^dz, 

a=- (d.Pdz + dtPdt + -dt] , (3=- ( d^Pdt + dtPdz + -dz 



t J ' ' 2 V t 
7 = i Qdt - dtPdt - d^Pdz^ , 6=^ (^jdz - dtPdz - d^Pdt 

(26) ^ ^ \ (d^Xdt + dtXdz - ^dz 
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where for the derivatives of P and A we get 



dtP{t,z) = -27rJi(27rt)cos27rz, 
d:,P{t,z) = -27rJo(27rt)sin27rz, 
dtX{t, z) = t {{dtP{t, z)f + (d,P{t, z))2) , 



(27) d,X{t, z) = 2t {dtP{t, z)) {d,P{t, z)) . 



5.2. Association with discrete forms. The next step is to choose an initial 
hypersurface. This has to obey our requirement that its topology is S^. That 
means we need a closed space-like curve in the two-dimensional orbit space 
For the test we used curves, whose i-coordinate is constant 



We subdivide this curve into n pieces, which are again of the form (|28p . but with 
A G [(i — l)/n,i/n],i = 1, . . . ,n. On each of these pieces we integrate the 1-forms, 
and on the nodes we evaluate the 0-forms to get initial values. 

Since light-rays in the polarized Gowdy solution ((25|l are the curves with = 
[dzj, the coordinate of the nodes in the computational mesh are easily obtained and 
it is thus straight forward to compare the numerical and the analytical solution 



This is done as follows. In every time slice there are n > 1 space-like edges. We 
take the maximum of the relative error of the invariant length^ of the edges in a 
time slice as a measure for the accuracy of the scheme. 

In this construction, since n is the number of initial edges, we need 2n time steps 
to simulate a crossing time {6t = 1). 



Now two scenarios in the polarised Gowdy solution are investigated, an expand- 
ing and a collapsing time evolution. For the expanding evolution we take an initial 
hypersurface of the form ([28]) with to = := 1 and zq — Zq := 0, and for the collaps- 
ing evolution these constants are to — '■— 9.8753205829098 and zq = Zq := —0.5. 

The hypersurface to = t^ is in the sense special that the Bessel function term 
Jo(27rt) vanishes and the Apples with Apples alliance suggests to use this particular 
choice. 

6.1. Expanding Evolution. In the expanding evolution we divide the initial curve 
into n edges of the form {z — X,t — 1, A G [{k — l)/n,k/n],k — 1, • . . , n}. For 
n = 50 and n = 100 the growth in time of the maximal relative error of the lengths 
of space-like edges is shown in figure [3] 

For scheme II we see that, although the absolute value of the lengths of the edges 
grows exponentially in this space-time, the relative error of these lengths grows only 
linearly and remains small even at t = 100. Furthermore in this scheme the error in 
the calculation with 100 initial edges is about 4 times smaller than the error with 
50 initial edges. This indicates quadratic convergence for decreasing edge lengths. 

"^The orbit space A4i can be identified with a surface that satisfies x = const, y = const. 
■^We define the relative error to be \lap/lex — 1|, where lap is the numerically calculated ap- 
proximation to the invariant length and lex is the length in the exact solution. 



(28) 




with A G [0, 1) and fixed to, zo- 



(see [28]). 



6. Concrete Examples and Results 
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Figure 3: Expanding evolution: Growth of the maximal relative error of the invari- 
ant lengths of the space-like edges. Left: Results of scheme I (logarithmic ordinate). 
Right: Results of scheme II (linear ordinate). 



For scheme I the error does not grow only linearly. We see that for three to four 
crossing times the scheme has a good numerical behaviour. The size of the error 
grows slowly, it is relatively small and becomes even smaller when more initial edges 
are used. But after four crossing times the error starts to grow exponentially, and 
the results of the calculations in different meshes become similar. Then, at i « 12 
the error in the simulation with 50 edges blows up rapidly. Similar behaviour can 
be observed for the simulation with 100 initial edges, but here the critical time is 
t « 20. 

In both schemes the growing of the error is superposed by a quasi-periodical 
fluctuation (for scheme II it is not visible in figure [3l because it is small, but it 
exists). The period of this fluctuation is St = 1/2, i.e. half as long as the crossing 
time. 

6.2. Collapsing Evolution. For the collapsing evolution we again divide the ini- 
tial curve into n edges of the form {z = X — 0.5, t = tg, A G [{k — l)/n, k/n], k = 
1, . . . , n}. For n = 50 and n = 100 the growth of the maximal relative error of the 
lengths of space- like edges is shown in figure ID 

Since at i = the polarised Gowdy space-time has a singularity (the cosmological 
singularity), we can only evolve as long as t remains positive. I.e. if we use for 
instance 50 initial edges, we can at most make 987 time-steps. 

A surprising result is that for both schemes the relative error does not grow at 
all. Instead it fluctuates quasi-periodically, but remains at similar sizes until the 
time approaches zero, i.e. near the singularity. 

Here the period of the fluctuation is St = 1/4, i.e. a quarter of a crossing time 
and only half as long as in the expanding time evolution. Further investigations 
reveal that for other initial hypersurfaces, i.e. to ^ tg, the error still fluctuates, but 
then the period of the fluctuation is 5t = 1/2, like in the expanding time evolution. 

Another point to mention is that for scheme II the error in the simulation with 
100 edges is about four times smaller than the error in the simulation with 50 
edges. This is not the case in scheme I where the errors of both simulations are of 
comparable size. 

6.3. Convergence of the Schemes. Finally we consider the convergence be- 
haviour of the numerical schemes for both examples. We saw that in the expanding 
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Figure 4: Collapsing evolution: Maximal relative error of the invariant lengths of 
the space-like edges. Left: Results of scheme I. Right: Results of scheme II. 
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Figure 5: Convergence of the schemes. Left: expanding evolution at t = 2. Right: 
collapsing evolution at t = — 4.9. 



time evolution scheme I is not stable, and that it behaves best for t < A. Therefore 
we investigate the errors sA t — 2. For the collapsing evolution we choose with 
t — Iq ~ 4.9 a time where the (fluctuating) error has neither a minimum nor a 
maximum. 

In figure [5] it is shown how both schemes behave when the typical length of the 
edges is decreased, i.e. when the number of initial edges becomes larger. To obtain 
the curves we calculated the time evolution with 10, 20, 40, 80 and 160 initial edges. 

We see that the results of scheme II converge quadratically to the analytical 
solution. The errors of scheme I on the other hand converge only linearly to zero 
in the expanding time evolution and do not converge at all in the collapsing time 
evolution. This is in agreement with the observations that we made in the time 
evolution simulations. 

6.4. Discussion. The first point to discuss is the quasi-periodical fluctuation of 
the error and its period of 6t — 1/2. This phenomenon points to the interpretation 
that, due to the periodic boundary conditions, the errors that sum up during the 
time evolution cancel out when a lightray intersects the other lightray that was sent 
from the same point in the other direction. 

In the collapsing time evolution we observed that the period of the fluctuation 
is only 5t — 1/4 when Iq = t§, but it is again 5t ~ 1/2 when other initial data are 
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chosen. This initial hypcrsurface is in the sense special that the Bessel function 
Jo(27rt) has a root there. Hence it is most likely that with these initial data other 
symmetries of the considered space-time cause new cancellations of errors. 

Concerning the quality of the numerical solutions we see that (at least for the 
considered examples) scheme II provides the best results. This scheme is accurate 
and stable, the relative error grows only linearly in time and remains small. Un- 
fortunately we cannot say much about its robustness (in the sense of the Apples 
with Apples alliance) , because we only considered the polarised Gowdy space-time. 
This scheme is the analogue of scheme III in [28], which also provided the best 
results there. The reason for this behaviour seems to be that the system (|8]) is 
geometrically preferred. 

The analogue of scheme I here is scheme I in [28]. For spherically symmetric 
space-times the results of that scheme were, compared to the results of scheme I 
here, better. Especially it was quadratically convergent. 

For scheme I here there are situations when it does not converge at all, and 
moreover it is not stable. We saw that for small times in the expanding evolution 
it behaves quite well, but after that there is a period when the error grows expo- 
nentially. Up to now it is not clear if the reason for this behaviour is rather the 
numerical scheme or the exponentially growing example. However, it becomes clear 
that the scheme is unstable when the period of exponential growth ends (at i « 12 
resp. t « 20). We thus conclude that scheme I is not practicable. 

Of course it is not too surprising that problems occur in scheme I, because the 
system ([7]) is not very well understood at the analytical level. In particular, we 
were not able to derive a minimal system for the gravitational wave space-times. 
Quite the contrary we found an overdetermined system, and although there are 
arguments that some equations are redundant, it was still necessary to choose 
discrete equations that we omitted without having justifying arguments for this 
particular choice. 

Comparing the previously presented numerical schemes with other codes one has 
to distinguish between the expanding and the collapsing case. The reason is that 
in the collapsing time evolution the coordinate representation (j25p is usually not 
used. Instead one takes a harmonic time coordinate r that becomes infinite at the 
singularity [25] . That means that even if the difference of t between the space-like 
slices is fixed, the difference of t becomes small. Yet, because of our method to 
construct the computational mesh we cannot control the slicing in this way. The 
results of other numerical methods are thus not really comparable. 

However, simulations of collapsing polarised Gowdy evolutions can be found e.g. 
in [12,25]. Furthermore [1] reports about results of the ADM and BSSN codes and 
the PITT group [23] presents results of the ABIGEL code. Further information 
about cosmological space-times with singularities can be found e.g. in [3]. 

Numerical simulations that use the expanding polarised Gowdy solution as a 
testbed can be found e.g. in [21]. Again the PITT group [23] presents results of 
the ABIGEL code. 

Summarising the hitherto existing work about discrete differential forms in nu- 
merical GR it seems that the most urgent task, in order to apply this method to 
physically more relevant situations, is to find statements about the underlying sys- 
tem of equations. Indeed it is not clear how to apply standard concepts of numerical 
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GR, like hyperbolicity, to the Cartan formulation. What is currently lacking is a 
better understanding of the notions of a discrete geometry. 
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Appendix A. Derivation oe the reduced system 

We are interested in space-times {A4 , g) which can be characterised by the fol- 
lowing conditions: 

(1) There is an effective, isometric action of on M. 

(2) The generators of this action, the two Killing vectors ^ and 77, are space-like, 
they are orthogonal and they commute. 

(3) The action is hypersurface-orthogonal, i.e., the 2-flats orthogonal to the 
span of the Killing vectors are integrable. 

From condition (1) it follows that the space-time is locally foliated by 2-dimensional 
orbits whose tangent spaces are spanned by the Killing vectors. Since these are 
orthogonal and space-like we can define two space-like unit vectors by 

(29) e = 6-^62, ry = eSe3, 

with functions / and g which are constant on the orbits. We complete these vec- 
tors to an orthonormal basis (60,61,62,63) with dual basis {9'^ ,9^ ,6^ ,6^). Since 
[62 , 63] = we have 

(30) d0°(62, 63) - 62(0^(63)) - 63(00(62)) - 0*'([e2, 63]) - 

and, similarly, for 6^. Hence, we get the expansions 

(31) d9°^Ae° + Be\ d6>i = C0" + D6>i 

for some 1-forms A, B, C and D. This can be expressed by the well known 
integrability conditions 

(32) de^9"9^ = 0, de^9"9^ = 0. 

While the choice of 62 and 63 is fixed by aligning them with the Killing vectors the 
other two basis vectors are not fixed. They can still rotate inside the plane which 
they generate with a rotation which may not only depend on the orbit but even the 
location on the orbit. We can fix the latter dependence by requiring that 60 and 
61 commute with both Killing vectors. This requirement is consistent because the 
Killing vectors commute. 

Then the four basis vectors are Lie transported along both Killing vectors and 
the same is true for the dual basis vectors, i.e., the equations 

(33) £^6" = = £^9" 

hold for a = 0, ... ,3. Using these equations and the expansions ([3T|) leads us after 
a short calculation to the fact that the 1-forms A to D are linear combinations of 
9^ and 9^ only, i.e., they vanish when restricted to an orbit. 
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Furthermore, condition (3) imphes that the distribution generated by eo and ei 
is integrable, i.e., that 

(34) [eo,ei] = aeo + l3ei. 

This imphes that there exist 2-dimensional submanifolds orthogonal to the orbits 
of the action which are dragged into each other by the group action. 
Expressed in terms of the dual basis, the integrability condition is 

(35) do'^e^e^ = 0, do^e^e^ = o, 

which leads to a similar expansion of the differentials in terms of the coframe as 
above. In an analogous calculation making use of the invariance of the frame (j33[) 
we find the following relations 

(36) de^ = dfe^, de^ = dge^. 

These relationships allow us to partly determine the connection forms from the first 
structure equation 

(37) d6>° + a;%6»'' = 0. 

We find that all the connection forms except for u; := uj^i are determined from the 
above equations 

(38) u^\^foe\ u;\^he\ Lj\^g^e\ u:\^gie\ c^2^ = 0, 

where /o = eo(/), etc. so that df = fo9^ + fi9^. So we find that the information 
about the extrinsic geometry of the orbits is contained in the two scalars / and 
g, while the geometry of the orbit space A4i which can be identified with the 
space orthogonal to the orbits is described by the coframe vectors (0*^, 6^) and the 
connection form w. 

With these expressions for the connection forms we get the Nester-Witten 2- 
forms 

Lo = {fi + gi)e^e\ L2 = -{^dg)e^ -u,e\ 

Li = ifo+go)e^e\ L3 = {*df)e^+u:e^ 



and the Sparling 3-forms 

So^{{fo+go)'^ + fidg + fo{*dg))e^e\ S2^u:dge\ 
Si = iih+gi)u + fodg + h{*dg))e^e\ S3 = -ujdfe\ 

Here, we have used the notation *q: for the 1-form aoO^ + ai9^ given the 1-form 
a = aod'^ + ai6^ . Inserting these expressions into the equations dLi = Si (and 
stripping off the common factor 9 9,9 and 9 , respectively) yields the four 
equations on A^i 

d(/i + gi) + ifi + .9i)(d/ + dg) = (/o + .go)'^ + /idg + /o {*dg), 
d(/o + go) + ih + .9o)(d/ + dg) = [h + gi)^ + fodg + h {^dg), 
^ ' d(*d.g) + dg(*d.g) + du; 0, 

d(*d/) + d/(*d/) + du; = 0. 

The first pair of equations are 1-form equations. They comprise four scalar equa- 
tions. However, they contain w which indicates that they depend on the chosen 
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gauge. We can extract two gauge invariant equations by multiplying with the basis 
1-forms and combining the equations linearly. This leads to the 2-form equations 

d2(/ + .g) + d(/ + 5)d(/ + .g) = 0, 

d(*d/) + d(*d5) + d(/ + .g)((*d/) + (*d.g)) = 0. 

The first of these is an identity while the second equation can be combined with the 
third and fourth equation of ((4T|) . So we have the following three 2-form equations 

dw = d/(*d5), 

d(*d5) + d.g(*d5) + d/(*d.g) - 0, 

d(*d/) + d/(*d/) + dg(*d/) = 0. 

Note, that we do not attach a meaning to the forms -kdf and -kdg as yet. These 
symbols are so far only abbreviations for the forms -kdf — /qB^ + and -kdg = 

Thus, we end up with the following set of differential forms which we regard as 
being defined on the frame bundle over A4i 

Zo = d(/o + go) + (/o + 5o)(d/ + dg) 

~ (A + g,)Lj - (/o.go + fi9i)e" - (Affo + fo9i)d\ 
Zi =d(A+gi) + (A+.9i)(d/ + d.g) 

- (/o + go)u: - (/o5o + fi9i)0' - {fm + h9i)0\ 
Z2 EEda;-d/(*d.g), 
Z3 = d{-kdg) + dg{-kdg) + d/(*dg), 
Z4 = d(*d/) + d/(*d/) + d,g(*d/), 

Z5 = de° + u}e^, 

Zg = d^i 

Z7 = d/ - /o0" - Ae\ 

Zg = dg - .go©" - .gif', 

Zg = (d/o - a;A)0" + (dA - ^k)0\ 
Zio = (dgo - t^gi)©" + (dgi - a;5o)0\ 
Zn = *d/ - /o0i - hO\ 
Z12 ^★dg-goe' - .gi0°, 

Zi3 = d(*d/) - (d/o - uh)e^ - (dA - ^/o)0°, 
Zi4 = d(*d.g) - (dgo - ^^91)^^ - (dgi - u>go)0° ■ 

Apart from the forms arising from the Witten- Sparling equations this system in- 
cludes forms which implement the expansion of d/, dg, -kdf and -kdg in terms of 
their components /o, A: 9o and gi as well as their exterior derivatives. 

Note, that there is also the first Bianchi identity which should be added. It reads 

(42) da;0° = = da;6»\ 

However, since it is an identity and, in addition, a 3-form equation which cannot 
be implemented in our 2-dimensional problem we will simply assume its general 
validity without listing it among the relevant forms. 
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A closer and somewhat tedious examination of the system reveals the following 
facts: 

dZo € (Zq, Zi, Z2, Z3, Z4, Z7, Zg, Zg, Zn, Z12, Z13) 
dZi G (Zq, Zi, Z2, Z3, Z4, Z7, Zg, Zio, Zn, Z12, Z14) 
dZ2 e (Z3) n (Z4) , dZ3 G (Z3) , dZ4 e (Z4) , 
dZ5e(Z6), dZ6e(Z5), 
dZ7e(Z9), dZ8e(Zio), dZ9G(Z5,Z6), dZm e (Z5,Z6) , 

dZnG(Zi3), dZi2e(Zi4), dZi3 e (Z5,Z6), dZu G (Z5, Zg) . 

The crucial forms are Z2, . . . , Zg. Z2, Z3, Z4 determine the connection and Z5 and 
Zg determine the metric in terms of the frame O'^ and 6^. The forms Z7, . . . , Z14 
are 'book keeping' forms in the sense that they express d/, dg, *d/ and -kdg in 
terms of the components /o,/i, ffo and gi. These are not really necessary if we are 
only interested in the geometry of A4i. In fact, in this case it is not even necessary 
to know about the functions / and g. 

We can extract several systems from these equations. The first system is ob- 
tained by ignoring that /o, ... ,171 are components of 1-forms with respect to the 
coframe and working simply with these four functions. The system we use consists 
of Zq, Zi, Z5, Zg, as well as Z3 + Z2 and Z4 + Z2. In addition, we include Zg, i.e., 
the exterior derivative of Z7. The corresponding equation for go and gi follows 
from Zq and Zi. This is the system ([T]). 

The second system is obtained by eliminating the components as completely as 
possible. We use generic 1-forms a, (3, 7 and d for the 1-forms d/, *d/, dg and 
★dg, respectively. Then the system consists of Zg^^ -I- Zi^^, linear combinations of 
Z2, Z3 and Z4 as well as Z5 and Zg. Finally we add the closedness of a and 7 and 
the algebraic relationships between a and f3 respectively 7 and S. This results in 
the system (|8]). 

Note, that it is enough to include the combination of Zq and Zi given above. 
This is due to the following consideration. Let / — (Zfc,2 < k < 14) be the ideal 
generated by all forms except for Zq and Zi. Then we have 

dZo -u;Zi +d(/ + .g)Zo = mod/, 
^'^'^^ dZi -u;Zo + d(/ + 5)Zi =0 mod/. 

By taking another exterior derivative we get the equations 

da; Zi = mod /, 
^"^^^ dwZoEEO mod/, 

since dl C /. This implies, that if all other equations hold, then we have duj Z^ = 
which forces Zi = unless du) = 0, a case in which we are not interested. 
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